Thalamocortical Connectivity Development along the Sensorimotor-Association Axis: Sensitivity Analyses

Read in Data for Developmental Characterization

Glasser parcel names

#Glasser regions with corresponding labels and tract names
glasser.regions <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360_regionlist.csv")
glasser.regions$tract <- paste0("thalamus_", glasser.regions$orig_parcelname) %>% gsub("_ROI", "_autotrack", .) %>% gsub("-", "_", .)

S-A axis

S.A.axis.cifti <- read_cifti("/cbica/projects/thalamocortical_development//Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo
S.A.axis <- data.frame(SA.axis = rank(S.A.axis.cifti$data), orig_parcelname = names(S.A.axis.cifti$Parcel))
S.A.axis <- merge(S.A.axis, glasser.regions, by = c("orig_parcelname"), sort = F)

Spatial permutation (spin) test parcel rotation matrix

perm.id.full <- readRDS("/cbica/projects/thalamocortical_development/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spins
spin.parcels <- glasser.regions

Dataset-specific tract lists

#generated by tract_measures/tractlists/thalamocortical_tractlists.R
tractlist.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractlist_primary.csv") 
tractlist.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractlist_primary.csv")

Sensitivity variable dfs

#Surface area
area.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/cortical_anatomy/individual/PNC/PNC_surfacearea_finalsample.csv")
area.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/cortical_anatomy/individual/HCPD/HCPD_surfacearea_finalsample.csv")

#SNR 
snr.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_SNR_finalsample.csv")
snr.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_SNR_finalsample.csv") 

#Overlap sensitivity
sensitivity.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_overlapsensitivity_finalsample.csv")
sensitivity.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_overlapsensitivity_finalsample.csv")

Sensitivity analysis results

setwd("/cbica/projects/thalamocortical_development/thalamocortical_results/sensitivity_results/") #gam output dir, from gam_functions/fit_sensitivityGams.R

#list all .csv files in the environment results directory
files <- list.files(getwd(), pattern = '.csv', ignore.case = T, full.names = F) 

for(i in 1:length(files)){
  filename <- files[i]
  Rname <- gsub('^.{12}|.{4}$', '', filename)
  x <- read.csv(filename, header = TRUE)
  assign(Rname, x) 
}
rm(x)

Association Between Tract-wise Sensitivity Variables and Main Analysis Maturational Timing

Main analysis thalamocortical development results

gameffects_FA_glasser_pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/development_gameffects_FA_glasser_pnc.csv")

gameffects_FA_glasser_hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/development_gameffects_FA_glasser_pnc.csv")

Sensitivity variables in thalamocortical connections

tract_average_measures <- function(input.df, average.measure, final.tractlist){
  input.df <- input.df %>% select(contains("autotrack")) #ensure we only have tract columns 
  tract.averages <- colMeans(input.df, na.rm = TRUE) %>% as.data.frame %>% set_names(sprintf("mean_%s", average.measure)) %>% mutate(tract = row.names(.))
  tract.averages <- tract.averages[tract.averages$tract %in% final.tractlist$tract,]
  tract.averages <- merge(tract.averages, S.A.axis, by = "tract", sort = F)
  return(tract.averages)
}

#surface area
regional.area.pnc <- tract_average_measures(input.df = area.glasser.pnc, average.measure = "area", final.tractlist = tractlist.pnc)
regional.area.hcpd <- tract_average_measures(input.df = area.glasser.hcpd, average.measure = "area", final.tractlist = tractlist.hcpd)

#connection SNR
tract.snr.pnc <- tract_average_measures(input.df = snr.glasser.pnc, average.measure = "SNR", final.tractlist = tractlist.pnc)
tract.snr.hcpd <- tract_average_measures(input.df = snr.glasser.hcpd, average.measure = "SNR", final.tractlist = tractlist.hcpd)

#atlas overlap sensitivity
tract.overlap.pnc <- tract_average_measures(input.df = sensitivity.glasser.pnc, average.measure = "sensitivity", final.tractlist = tractlist.pnc)
tract.overlap.hcpd <- tract_average_measures(input.df = sensitivity.glasser.hcpd, average.measure = "sensitivity", final.tractlist = tractlist.hcpd)

#Combine development results and sensitivity measures
df.list.pnc <- list(gameffects_FA_glasser_pnc, regional.area.pnc, tract.snr.pnc, tract.overlap.pnc) #dfs to merge
gameffects_FA_glasser_pnc <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list.pnc) 
gameffects_FA_glasser_pnc <- left_join(spin.parcels, gameffects_FA_glasser_pnc, by = join_by(orig_parcelname, label, tract))

df.list.hcpd <- list(gameffects_FA_glasser_hcpd, regional.area.hcpd, tract.snr.hcpd, tract.overlap.hcpd) #dfs to merge
gameffects_FA_glasser_hcpd <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list.hcpd) 
gameffects_FA_glasser_hcpd <- left_join(spin.parcels, gameffects_FA_glasser_hcpd, by = join_by(orig_parcelname, label, tract))

Associations

surface area

cor.test(gameffects_FA_glasser_pnc$smooth.increase.offset, gameffects_FA_glasser_pnc$mean_area, method = c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_pnc$smooth.increase.offset and gameffects_FA_glasser_pnc$mean_area
## S = 1273292, p-value = 0.5268
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.04500699

perm.sphere.p(x = gameffects_FA_glasser_pnc$smooth.increase.offset, y = gameffects_FA_glasser_pnc$mean_area, perm.id = perm.id.full, corr.type = "spearman")

## [1] 0.35995

cor.test(gameffects_FA_glasser_hcpd$smooth.increase.offset, gameffects_FA_glasser_hcpd$mean_area, method = c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_hcpd$smooth.increase.offset and gameffects_FA_glasser_hcpd$mean_area
## S = 1238635, p-value = 0.5516
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.04256293

perm.sphere.p(x = gameffects_FA_glasser_hcpd$smooth.increase.offset, y = gameffects_FA_glasser_hcpd$mean_area, perm.id = perm.id.full, corr.type = "spearman")

## [1] 0.37435

connection SNR

cor.test(gameffects_FA_glasser_pnc$smooth.increase.offset, gameffects_FA_glasser_pnc$mean_SNR, method = c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_pnc$smooth.increase.offset and gameffects_FA_glasser_pnc$mean_SNR
## S = 1365701, p-value = 0.7327
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.0243017

perm.sphere.p(x = gameffects_FA_glasser_pnc$smooth.increase.offset, y = gameffects_FA_glasser_pnc$mean_SNR, perm.id = perm.id.full, corr.type = "spearman")

## [1] 0.44055

cor.test(gameffects_FA_glasser_hcpd$smooth.increase.offset, gameffects_FA_glasser_hcpd$mean_SNR, method = c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_hcpd$smooth.increase.offset and gameffects_FA_glasser_hcpd$mean_SNR
## S = 1063898, p-value = 0.01229
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1776307

perm.sphere.p(x = gameffects_FA_glasser_hcpd$smooth.increase.offset, y = gameffects_FA_glasser_hcpd$mean_SNR, perm.id = perm.id.full, corr.type = "spearman")

## [1] 0.1184

atlas overlap sensitivity

cor.test(gameffects_FA_glasser_pnc$smooth.increase.offset, gameffects_FA_glasser_pnc$mean_sensitivity, method = c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_pnc$smooth.increase.offset and gameffects_FA_glasser_pnc$mean_sensitivity
## S = 1135068, p-value = 0.03563
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1486775

perm.sphere.p(x = gameffects_FA_glasser_pnc$smooth.increase.offset, y = gameffects_FA_glasser_pnc$mean_sensitivity, perm.id = perm.id.full, corr.type = "spearman")

## [1] 0.14175

cor.test(gameffects_FA_glasser_hcpd$smooth.increase.offset, gameffects_FA_glasser_hcpd$mean_sensitivity, method = c("spearman"))

## 
##  Spearman's rank correlation rho
## 
## data:  gameffects_FA_glasser_hcpd$smooth.increase.offset and gameffects_FA_glasser_hcpd$mean_sensitivity
## S = 1065814, p-value = 0.01305
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.1761501

perm.sphere.p(x = gameffects_FA_glasser_hcpd$smooth.increase.offset, y = gameffects_FA_glasser_hcpd$mean_sensitivity, perm.id = perm.id.full, corr.type = "spearman")

## [1] 0.10545

Developmental Sensitivity Analyses

Number of significant developmental effects

Function to calculate FDR-corrected significant age effects

sig.effects <- function(dataset, sensitivity_measure){
 
   ageeffects.df <- get(sprintf("gameffects_FA_%s_%s", sensitivity_measure, dataset)) 
  ageeffects.df$significant <- p.adjust(ageeffects.df$Anova.smooth.pvalue, method = c("fdr")) < 0.05 
  
  sigeffect.totaln <- sum(ageeffects.df$significant) 
  sigeffect.percent <- round(sigeffect.totaln/length(ageeffects.df$significant)*100, 2) 
  
  print(sprintf("In %s, %s thalamocortical connections (%s percent) show significant age-related changes when controlling for %s", dataset, sigeffect.totaln, sigeffect.percent, sensitivity_measure))
}

Sensitivity analysis: surface area

sig.effects(dataset = "pnc", sensitivity_measure = "area")

## [1] "In pnc, 202 thalamocortical connections (90.58 percent) show significant age-related changes when controlling for area"

sig.effects(dataset = "hcpd", sensitivity_measure = "area")

## [1] "In hcpd, 180 thalamocortical connections (78.26 percent) show significant age-related changes when controlling for area"

Sensitivity analysis: connection SNR

sig.effects(dataset = "pnc", sensitivity_measure = "SNR")

## [1] "In pnc, 201 thalamocortical connections (90.13 percent) show significant age-related changes when controlling for SNR"

sig.effects(dataset = "hcpd", sensitivity_measure = "SNR")

## [1] "In hcpd, 179 thalamocortical connections (77.83 percent) show significant age-related changes when controlling for SNR"

Sensitivity analysis: atlas overlap sensitivity (true positives)

sig.effects(dataset = "pnc", sensitivity_measure = "sensitivity")

## [1] "In pnc, 203 thalamocortical connections (91.03 percent) show significant age-related changes when controlling for sensitivity"

sig.effects(dataset = "hcpd", sensitivity_measure = "sensitivity")

## [1] "In hcpd, 177 thalamocortical connections (76.96 percent) show significant age-related changes when controlling for sensitivity"

Cortex-wide thalamic connectivity developmental profiles

Function to plot tract-wise GAM smooths

developmental.smooths <- function(dataset, sensitivity_measure){
 
  #get data for plotting 
  ageeffects.df <- get(sprintf("gameffects_FA_%s_%s", sensitivity_measure, dataset)) 
  smooths.df <- get(sprintf("smoothcentered_FA_%s_%s", sensitivity_measure, dataset)) 
  
  smooths.df <- merge((smooths.df %>% filter(grepl("thalamus_R", tract))), (ageeffects.df %>% select(tract, GAM.smooth.partialR2)), by = "tract")
  
  #match plotting params used in the main manuscript figures
  if(dataset == "pnc"){
    limit.low = 0
    limit.high = 0.1
  }
  if(dataset == "hcpd"){
    limit.low = -0.02
    limit.high = 0.075
  }
  
  #plot tract-wise developmental trajectories
  smooths.plot <- ggplot(smooths.df, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) +
    geom_line(linewidth = .3, alpha = .8) + 
    scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(limit.low, limit.high), oob = squish) +
    theme_minimal() +
    labs(x="\nAge", y="Connection FA\n") +
    scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
    theme(
    legend.position = "none", 
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text = element_text(size = 5, family = "Arial", color = c("black")),
    axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
    axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
    axis.line = element_line(linewidth = .2), 
    axis.ticks = element_line(linewidth = .2))
  
  ggsave(filename = sprintf("/cbica/projects/thalamocortical_development/figures/images/Figure_sensitivity_supplement/agesmooths_%s_%s.pdf", sensitivity_measure, dataset), device = "pdf", dpi = 500, width = 1.85, height = 1.65)

  return(smooths.plot)
}

Sensitivity analysis: surface area

developmental.smooths(dataset = "pnc", sensitivity_measure = "area")

developmental.smooths(dataset = "hcpd", sensitivity_measure = "area")

Sensitivity analysis: connection SNR

developmental.smooths(dataset = "pnc", sensitivity_measure = "SNR")

developmental.smooths(dataset = "hcpd", sensitivity_measure = "SNR")

Sensitivity analysis: atlas overlap sensitivity (true positives)

developmental.smooths(dataset = "pnc", sensitivity_measure = "sensitivity")

developmental.smooths(dataset = "hcpd", sensitivity_measure = "sensitivity")

Thalamocortical Connectivity Maturation is Organized by the Sensorimotor-Association Axis

Function to compute age of maturation correlation with the S-A axis and its spin-based significance

SA.alignment.statistics <- function(dataset, sensitivity_measure){
 
   ageeffects.df <- get(sprintf("gameffects_FA_%s_%s", sensitivity_measure, dataset)) 
   ageeffects.df <- merge(ageeffects.df, (S.A.axis %>% select(tract, SA.axis)), by = "tract", sort = F)
   
   spin.data <- left_join(spin.parcels, ageeffects.df, by = c("tract")) 
   corr.value <- cor.test(spin.data$SA.axis, spin.data$smooth.increase.offset, method = c("spearman"), exact = F)$estimate
   pspin.value <- perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman")
   
   results <- cbind(corr.value, pspin.value)

  return(results)
}

Function to plot thalamocortical connection age of maturation along the S-A axis

SA.alignment.plot <- function(dataset, sensitivity_measure){
 
   ageeffects.df <- get(sprintf("gameffects_FA_%s_%s", sensitivity_measure, dataset)) 
   ageeffects.df <- merge(ageeffects.df, (S.A.axis %>% select(tract, SA.axis)), by = "tract", sort = F)

   SA.plot <- ggplot(ageeffects.df, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
    geom_point(size = 0.5) +
    geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) + 
    scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
    theme_classic() +
    scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
    labs(x="\nS-A axis", y="Age of maturation (years)\n") +
    theme(
    legend.position = "none", 
    axis.text = element_text(size = 5, family = "Arial", color = c("black")),
    axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
    axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
    axis.line = element_line(linewidth = .2), 
    axis.ticks = element_line(linewidth = .2)) 

    ggsave(filename = sprintf("/cbica/projects/thalamocortical_development/figures/images/Figure_sensitivity_supplement/SAaxis_alignment_%s_%s.pdf", sensitivity_measure, dataset), device = "pdf", dpi = 500, width = 1.85, height = 1.65)

  return(SA.plot)
}

Sensitivity analysis: surface area

SA.alignment.statistics(dataset = "pnc", sensitivity_measure = "area")

##     corr.value pspin.value
## rho  0.5114256       6e-04

SA.alignment.plot(dataset = "pnc", sensitivity_measure = "area")

SA.alignment.statistics(dataset = "hcpd", sensitivity_measure = "area")

##     corr.value pspin.value
## rho  0.4913446     0.00185

SA.alignment.plot(dataset = "hcpd", sensitivity_measure = "area")

Sensitivity analysis: connection SNR

SA.alignment.statistics(dataset = "pnc", sensitivity_measure = "SNR")

##     corr.value pspin.value
## rho  0.5046296     0.00085

SA.alignment.plot(dataset = "pnc", sensitivity_measure = "SNR")

SA.alignment.statistics(dataset = "hcpd", sensitivity_measure = "SNR")

##     corr.value pspin.value
## rho  0.5129782     0.00215

SA.alignment.plot(dataset = "hcpd", sensitivity_measure = "SNR")

Sensitivity analysis: atlas overlap sensitivity (true positives)

SA.alignment.statistics(dataset = "pnc", sensitivity_measure = "sensitivity")

##     corr.value pspin.value
## rho  0.5293327     0.00055

SA.alignment.plot(dataset = "pnc", sensitivity_measure = "sensitivity")

SA.alignment.statistics(dataset = "hcpd", sensitivity_measure = "sensitivity")

##     corr.value pspin.value
## rho  0.5601407      0.0011

SA.alignment.plot(dataset = "hcpd", sensitivity_measure = "sensitivity")